The seismic attenuation effect needs to be taken into account for characterization of rock properties and proper amplitude-variations-with-offset (AVO) analysis. In migration, seismic attenuation information is needed to compensate for the absorption effect to enhance migration image quality. Therefore, estimation of seismic attenuation is essential for reservoir detection and monitoring.
Seismic attenuation can be quantitatively described by the quality factor Q. A simple assumption is that the seismic attenuation is frequency dependent but the quality factor Q is frequency independent. This assumption is valid in the frequency range of exploration geophysics applications. Q tomography is an approach for Q estimation. This approach attempts to reconstruct subsurface 2D or 3D Q models from seismic data. Generally, Q tomography algorithms are classified into two main categories. One category is ray-based tomography (Quan and Harris, 1997; Plessix, 2006; Rossi et al., 2007). The other category is wave-equation-based tomography (Liao and McMechan, 1996; Hicks and Pratt, 2001, Pratt et al., 2003; Watanabe et al., 2004; Gao et al., 2005). Wave-equation-based tomography is physically more accurate but computationally expensive and not practical for 3D cases. The present invention belongs to the category of ray-based Q tomography.
One major problem with Q tomography is how to establish the link between Q models and seismic data with minimum approximations and with maximum flexibility. A widely used approach is based on the relationship between Q and seismic amplitude decay. Another approach uses the seismic centroid frequency downshift to estimate the quality factor Q. This latter approach is believed to be more robust because this approach is independent of the geometrical spreading effect and reflection/transmission loss. However, the conventional centroid frequency shift method can use only Gaussian, boxcar, or a triangular function to fit the source amplitude spectrum, which introduces significant error because, in most cases, the source spectrum cannot be approximated by these functions. The present invention includes a frequency weighted exponential function, which is designed to fit various asymmetric source amplitude spectra to improve the accuracy of Q tomography by greatly reducing the source amplitude spectrum fitting error.
In most existing Q tomography algorithms, the optimization part is based on unconstrained optimization methods or based on simple nonnegative constrained optimization methods. As a result, these Q tomography algorithms take a large amount of computation time or yield many artifacts and unrealistic Q models (e.g., negative Q values or extremely low Q values), especially when the seismic data are contaminated by noise. The present invention, including an efficient optimization algorithm with box constraints, is able to improve the quality and the reliability of the reconstructed Q models. A more detailed discussion of the prior art follows next.
Seismic attenuation tomography (Q tomography) has been investigated for many years and much progress has been made. The two main components of a ray-based Q tomography algorithm are 1) a simple but accurate relationship between seismic data and Q models for constructing the mathematical model for Q tomography; 2) a reliable and robust optimization algorithm for solving this mathematical problem. Many techniques were developed or proposed for building these two components. These techniques are discussed below.
Establishing Link between Seismic Data and Q Models
The most simple and straightforward method to estimate Q is the spectral ratio method (Spencer et al., 1982; Tonn, 1991), where the logarithm of the spectra ratio between two seismic waveforms is calculated as a function of frequency, and this function is approximated by a linear function of frequency, whose slope is treated as the accumulated seismic attenuation and is eventually related to the Q values along the wave propagation path. Ideally, this method removes the effect of geometrical spreading and reflection/transmission loss with the assumption that these effects are frequency independent. In practical applications, this method is relatively unreliable due to wavelet overlapping, uncertainty in linear fitting, and many other factors.
Rickett (2006) proposed a tomographic extension of the spectral ratio method with the aid of time-frequency analysis technique. This approach was claimed to be insensitive to absolute scaling and was applied in an application of Q profile estimation using a vertical seismic profile (VSP). In this approach, the log-amplitude scalars describing the frequency independent amplitude variation are included in the unknowns, which substantially increases the number of unknowns and reduces the efficiency of the algorithm. Furthermore, in 2D/3D Q tomography using surface seismic reflection data, the log-amplitude scalars are not only a function of position, but also a function of ray, which severely complicates the procedure.
Based on the fact that variation of seismic wavelet rise-time is linearly related to the 1/Q profile along the propagation path, Wu and Lees (1996) reported a seismic attenuation tomography method using the rise-time in earthquake seismology. Unfortunately, this method is impractical in exploration geophysics because the wavelets are inevitably contaminated by noise, scattering effect, overlapping, etc.
It was pointed out that the shape of the seismic wavelet amplitude spectrum is almost exclusively affected by the quality factor Q, and a peak frequency variation method was developed for Q estimation (Zhang, 2008). This method is attractive but, in practice, there are difficulties in accurately measuring peak frequency variation. Moreover, because only the information at an individual frequency is used, the uncertainty of Q estimation can be large.
A more robust method was introduced by Quan and Harris (1997), where the information over the whole frequency band of seismic waveforms is used to calculate the centroid frequency downshift and then relate the centroid frequency shift to the Q profile along the raypath by a simple closed form formulation. This method is intrinsically immune to geometrical spreading and reflection/transmission loss. The limitation of this method is that the source amplitude spectrum has to be a Gaussian, boxcar, or triangular function. It is well known that the seismic amplitude spectrum is never a boxcar or triangular function. Also, it is usually asymmetric and can be very different from a Gaussian function. If this asymmetric amplitude spectrum is approximated by a Gaussian function, significant errors will be introduced in the reconstruction of Q models. Therefore, if there is a function that can be used to fit various seismic frequency spectra more accurately without losing the simple nature of the relation between the centroid by frequency of the recorded seismic data and the Q profile along the raypath, this robust method can be more accurate in practical Q tomography applications.
Constrained Optimization Algorithms for Q Tomography
When the relation between seismic data and Q models is established, the ray-based Q tomography problem can be described by a linear optimization problem. In most existing Q tomography algorithms, this linear optimization problem is solved iteratively by using Krylov subspace methods, such as the conjugate gradient method and the LSQR method without applying any constraints (Quan and Harris, 1997; Plessix 2006, Rossi et al, 2007). These algorithms work well provided that the seismic data have high signal-to-noise ratio (SNR). However, seismic data are never clean; in handling real field data, these unconstrained optimization algorithms always result in some negative Q values or extremely small positive Q values, which are physically unreal. Furthermore, under some circumstances, a priori information of the range of Q values is known. In these cases, this information needs to be included in the tomography algorithm through a box constraint to provide more reliable Q tomography results.
Rickett (2006) developed a Q estimation algorithm with a constraint. But his algorithm adopts a nonnegative constraint instead of a box constraint, which means, negative Q values are eliminated but extremely low positive Q values can still exist. Rickett (2006) reported two methods to apply the nonnegative constraint. The first method is a nonlinear transformation method. In this nonlinear transformation method, the variable Q is replaced by ey and y is solved for instead of 1/Q. By doing this, Q is forced to be positive but the whole system can be very nonlinear. To achieve this goal, the resulting optimization system is solved with a Gauss-Newton approach, which can be very expensive. Another disadvantage of this method is that, during optimization, when the Q values are very large, the gradient-based optimization algorithm will be stagnant or converge very slowly, i.e., the Q values will stay there and no longer change. In the worst case, if the y values in the starting model are infinite, then the gradient of the cost function is 0 and the optimization algorithm does not perform searching. The second method of applying the nonnegative constraint reported by Rickett (2006) is the enforcement of the monotonically increasing property of attenuation by a smoothing technique. This method works effectively for Q estimation using VSP data but may fail in 2D Q tomography using surface reflection seismic data.
The present inventors know of no existing Q tomography algorithm with box constraints to enforce the estimated Q values within the ranges specified by the upper boundaries and the lower boundaries. However, optimization algorithms with box constraints are employed in some other geophysical applications, such as velocity tomography. For example, Delbos et al. (2006) developed a seismic reflection tomography algorithm with box constraints. In their algorithm, the constrained optimization problem is solved with a Gauss-Newton augmented Lagrangian approach and the associated Lagrange problem, another constrained optimization problem, is solved by a combination of the gradient projection method, the active-set method, and the conjugate gradient method. The active-set method they use is conventional, and it is inefficient because the algorithm updates the active set, one constraint at a time (Bierlaire et al., 1991). When the number of box constraints is huge, the convergence rate of the algorithm can be very slow.
In the present invention, a recent development in the mathematical field, which may be referred to as the multi-index active-set method (Morigi et al., 2007), is employed to perform the Q tomography with box constraints, which significantly improves the performance of the Q tomography algorithm in terms of Q reconstruction quality and algorithm efficiency compared with the unconstrained Q tomography algorithm and the constrained algorithm using the conventional active-set method.